HPC_loop <- function(s){
  
  set.seed((2000  + s))
  design_num <- s %% nrow(design)
  if (design_num == 0){design_num <- nrow(design)}
  
  # Read in from design
  ss_use <- design[design_num, 1]
  vio_use <- design[design_num, 2]
  div_num <- design[design_num, 3]
  choice_divergence <- choice_divergences[div_num]
  outcome_divergence <- outcome_divergences[div_num]
  
  output.path <- file.path(paste0(use_path, s, ".rds"))
  
  bounds.sim <- tryCatch({
    
    # ###############
    # Original PICA 
    # ###############
    d.sim_pica1 <- make_sim_data(ss_use, choice_divergence, outcome_divergence)
    
    
    effects.focus <- rbind(c(1, 3, -1), c(2, 3, -2))
    colnames(effects.focus) <- c("treat", "ctrl", "C")
    
    bounds_1 <- bounds.frechet_new_estimand(data = d.sim_pica1$data, J = 3,
                                            effects.mat = effects.mat,
                                            effects.focus = effects.focus,
                                            posterior = TRUE,  
                                            n_dir_MC = n_dir_MC,
                                            n_norm_MC = n_norm_MC,
                                            hpd = TRUE,
                                            alpha = .95,
                                            sensitivity = FALSE,
                                            rhos = Inf
    )
    
    bounds.est_1  <- bounds_1
    bounds.post_1 <- bounds_1$posterior
    bounds.out_1 <- cbind(bounds_1$effects, bounds_1$posterior$effects.ci[, -c(1,2,3)])
    bounds.out_1 <- rbind(bounds.out_1, bounds_1$effects.new_tab)
    
    bounds.boot_1 <- NULL
    
    # ##############
    # PICA-2
    # ##############
    d.sim <- make_sim_data_placebo(n = ss_use, choice_divergence, outcome_divergence, 
                                   vio = vio_use)
    
    # Analyze as the PICA-2 
    out_2 <- pica2(Y = "Y", D = "D", A = "A", C = "C", 
                   placebo = "pl", placebo_level = 3, data = d.sim$data)
    
    out_2_tab0 <- rbind(out_2$ATE, out_2$ACTE_view, out_2$ACTE_non_view)
    C     <- c(0, 0, 1, 2, -1, -2)
    treat <- c(1, 2, 1, 2, 1, 2)
    ctrl  <- c(3, 3, 3, 3, 3, 3)
    out_2_tab <- cbind(C, treat, ctrl, out_2_tab0)
    
    test_2 <- test_pica2(Y = "Y", D = "D", A = "A", C = "C", 
                         placebo = "pl", placebo_level = 3, data = d.sim$data)
    
    test_2_boot <- test_pica2_boot(Y = "Y", D = "D", A = "A", C = "C", 
                                   placebo = "pl", placebo_level = 3, data = d.sim$data,
                                   num_boot = 100)
    
    test_2_res <- c(test_2["diff"], test_2["diff_SE"], test_2["p-value"],
                    test_2_boot["diff"], test_2_boot["diff_SE"], test_2_boot["p-value"])
    
    # Analyze as the usual PICA
    out_2_1 <- pica1_after_pica2(Y = "Y", D = "D", A = "A", C = "C", S = "S",
                                 placebo = "pl", placebo_level = 3, data = d.sim$data,
                                 posterior = TRUE,
                                 n_dir_MC = n_dir_MC,
                                 n_norm_MC = n_norm_MC,
                                 hpd = TRUE,
                                 alpha = .95,
                                 sensitivity = FALSE,
                                 rhos = Inf)
    
    out_2_1_boot <- NULL
    
    list(design_num = design_num, 
         bounds.out_1 = bounds.out_1,
         bounds.est_1 = bounds.est_1,
         bounds.post_1 = bounds.post_1,
         out_2_tab = out_2_tab,
         test_pica2 = test_2_res,
         out_2_1_tab = out_2_1,
         out_2_1_boot_tab = out_2_1_boot
    )
    
  },
  error = function(e){
    message('caught error in draw ', s, ' (n=', ss_use, '): ', as.character(e))
  })
  
  saveRDS(bounds.sim, output.path)
}